rm(list = ls())
warnings()
library(fOptions)

########################################################
#I. MINISUM FUNCTION----THEOREM 1---INDEPENDENT PROGRAM#
########################################################
f_h<-function(sigma,P){

print('inside f_h')
#computes function h(sigma)=C_BS(X1, sigma)-C_BS(X2,sigma) 
#required for the theorem 1 and function f_minisum

Call_mkt_1=unlist(P[1])
Call_mkt_2=unlist(P[2])
maturity_dt=unlist(P[3])
X1=unlist(P[4])
X2=unlist(P[5])
S=unlist(P[6])
risk_free=unlist(P[7])
dividend=unlist(P[8]) # cost of carry  or dividends,take them zero in most cases

# in the function help it is explained that to get BS formula we have to set b=risk_free
# to see the content of the function enter GBSOption

Call_BS_1=GBSOption(TypeFlag = "c", S = S, X = X1, Time =maturity_dt, r = risk_free,  b =risk_free, sigma = sigma)
Call_BS_2=GBSOption(TypeFlag = "c", S = S, X = X2, Time =maturity_dt, r = risk_free,  b =risk_free, sigma = sigma)

h_out=Call_BS_1@price-Call_BS_2@price
cat('h_out=',h_out,'\n')

return(h_out)
}#end of f_h

#################################################################

f_f<-function(sigma, P){
#computes function f(sigma) defined in the paper, its minimum gives minisum

Call_mkt_1=unlist(P[1])
Call_mkt_2=unlist(P[2])
maturity_dt=unlist(P[3])
X1=unlist(P[4])
X2=unlist(P[5])
S=unlist(P[6])
risk_free=unlist(P[7])
dividend=unlist(P[8]) # cost of carry  or dividends,take them zero in most cases

# in the function help it is explained that to get BS formula we have to set b=risk_free
# to see the content of the function enter GBSOption

Call_BS_1=GBSOption(TypeFlag = "c", S = S, X = X1, Time =maturity_dt, r = risk_free,  b =risk_free, sigma = sigma)
Call_BS_2=GBSOption(TypeFlag = "c", S = S, X = X2, Time =maturity_dt, r = risk_free,  b =risk_free, sigma = sigma)

f=abs(Call_mkt_1 - Call_BS_1@price)+abs(Call_mkt_2-Call_BS_2@price)
return(f)
} # end of f_function

##########################################

f_minisum<-function(ls_ParamOptions ){
# computes both the minisum function , i.e. min of the f_function defined above, and
# the arguments of f_function  at which the function reaches min, i.e. sigma_star
#requires functions: f_f and f_h
#ls_ParamOptions=list(Call_mkt_1,Call_mkt_2, maturity_dt, X1,X2,S,risk_free, dividend) 


print('inside f_minisum')

Call_mkt_1=unlist(ls_ParamOptions[1])
Call_mkt_2=unlist(ls_ParamOptions[2])
maturity_dt=unlist(ls_ParamOptions[3])
X1=unlist(ls_ParamOptions[4])
X2=unlist(ls_ParamOptions[5])
S=unlist(ls_ParamOptions[6])
risk_free=unlist(ls_ParamOptions[7])
dividend=unlist(ls_ParamOptions[8])

length_data=length(Call_mkt_1)

sigma_star<-vector("numeric",length = length_data)
flag<-vector("numeric", length = length_data)
Minisum<-vector("numeric", length = length_data)
fi<-vector("numeric", length = length_data)


for (i in 1:length_data) {


 sigma_star[i]=0
 Minisum[i]=0
 flag[i]=0 # to check that one of the conditions below was satisfied, in that case flag[i] becomes 1
 
 ls_ParamOptions_i=list(Call_mkt_1[i],Call_mkt_2[i], maturity_dt[i], X1[i],X2[i],S[i],risk_free[i], dividend[i]) 

 fi[i]=(log(X1[i]*X2[i]/(S[i]*S[i]))/maturity_dt[i])-2*risk_free[i]

 sigma_imp1=GBSVolatility(price=Call_mkt_1[i], TypeFlag="c", S=S[i], X=X1[i], Time=maturity_dt[i], r=risk_free[i], b=risk_free[i], tol=0.0001, maxiter=100)
 sigma_imp2=GBSVolatility(price=Call_mkt_2[i], TypeFlag="c", S=S[i], X=X2[i], Time=maturity_dt[i], r=risk_free[i], b=risk_free[i], tol=0.0001, maxiter=100)

cat( ' sigma_imp1=',sigma_imp1,'\n')
cat(' sigma_imp2=', sigma_imp2,'\n')

 f1=f_f(sigma_imp1,ls_ParamOptions_i)
 f2=f_f(sigma_imp2,ls_ParamOptions_i)

      #Compare sigma_imp1, sigma_imp2 and fi_sqrt to determine value of minisum at sigma_star
	if(fi[i]<=0){
	 Minisum[i]=f2
	 sigma_star[i]=sigma_imp2
	 flag[i]=1.1+flag[i]

	 next       # forces to start a new iteration of the loop
	} # end if(fi[i]<=0)

	if(fi[i]>0 ) fi_sqrt=sqrt(fi[i])

	if( (fi[i]>0) && (fi_sqrt<= min(sigma_imp1, sigma_imp2 )) ){
	 Minisum[i]=f2
	 sigma_star[i]=sigma_imp2
	 flag[i]=1.2+flag[i]
	 next
	}# end if(fi[i]>0 && fi_sqrt<= min(sigma_imp1, sigma_imp2 ))

	if( (fi[i]>0) && (sigma_imp1 < fi_sqrt) && (fi_sqrt < sigma_imp2)){
	 Minisum[i]=min(f1,f2)
	 #cat('Inside function f_minisum: now will be computing sigma_star=argmin( h(), h() \n')
	 #cat('sigma_imp1=',sigma_imp1  ,  'sigma_imp2=',sigma_imp2,'\n')
                  if( f_h(sigma_imp1,ls_ParamOptions_i)<f_h(sigma_imp2,ls_ParamOptions_i) ) {
                                                             sigma_star[i]= sigma_imp1} else {
                                                                    sigma_star[i]= sigma_imp2} #end if( f_h(sigma_imp1,ls_ParamOptions_i)<....

cat('Inside function f_minisum' ,'\n')
	 flag[i]=1.3+flag[i]
	 next
	}# end if( (fi>0) && (sigma_imp1 < fi_sqrt) && (fi_sqrt < sigma_imp2)){

	if((fi[i]>0) && (sigma_imp2 < fi_sqrt) && (fi_sqrt < sigma_imp1) ){
	 Minisum[i]=f_f(fi_sqrt,ls_ParamOptions_i)
	 sigma_star[i]=fi_sqrt
	 flag[i]=1.4+flag[i]
	 next
	}#end if((fi[i]>0) && (sigma_imp2 < fi_sqrt) && (fi_sqrt < sigma_imp1) )

	if(  (fi[i]>0) && (fi_sqrt >max(sigma_imp1,sigma_imp2) )  ){
	 Minisum[i]=f1
	 sigma_star[i]=sigma_imp1
	 flag[i]=1.5+flag[i] 
	 next
	}# end if(((fi[i]>0) && (fi_sqrt >max(sigma_imp1,sigma_imp2) )  )


  stop('Error in the function f_minisum: this line must not be reached ')

  }# end for (i in 1:length(Call_mkt_1))

  ls_out=list(Call_mkt_1, Call_mkt_2, maturity_dt, X1,X2,S,risk_free, dividend, Minisum,sigma_star, flag, fi)


return(ls_out)
}# end f_minisum

#######################################################
f_ReadData<-function(file_location_name){
# This function reads data, which location and name are given in the string file_location_name
# put data into several vectors and return a list of those vectors

print('f_ReadData: before data reading')
d<-read.table(file_location_name, header=TRUE, sep="," , na.string="." )
print('f_ReadData: after data reading')


length_data=length(d$Call_mkt_1)
print(length_data)


v_Call_mkt_1<-vector("numeric",length = length_data)
v_Call_mkt_2<-vector("numeric",length = length_data)
v_X1<-vector("numeric",length = length_data)
v_X2<-vector("numeric",length = length_data)
v_maturity_dt<-vector("numeric",length = length_data)
v_S<-vector("numeric",length = length_data)
v_risk_free<-vector("numeric",length = length_data)
v_dividend<-vector("numeric",length = length_data)


for (i in 1:length_data) { 

# Require X1<X2
	if(d$X1[i]<=d$X2[i]){
		v_Call_mkt_1[i]=d$Call_mkt_1[i]
		v_Call_mkt_2[i]=d$Call_mkt_2[i]
		v_X1[i]=d$X1[i]
		v_X2[i]=d$X2[i]
	}  else { # in case X1>X2 we switch the options

		v_Call_mkt_2[i]=d$Call_mkt_1[i]
		v_Call_mkt_1[i]=d$Call_mkt_2[i]
		v_X2[i]=d$X1[i]
		v_X1[i]=d$X2[i]
	}# end if(d$X1[i]<=d$X2[i])

v_maturity_dt[i]=d$maturity_dt[i]
v_S[i]=d$S[i]
v_risk_free[i]=d$risk_free[i]
v_dividend[i]=d$dividend[i]

ls_ParamOptions=list(v_Call_mkt_1, v_Call_mkt_2, v_maturity_dt, v_X1,v_X2,v_S,v_risk_free, v_dividend) 
} #end for (i in 1:length_data)

print('f_ReadData: ls_ParamOptions=')
print(ls_ParamOptions)
return(ls_ParamOptions)

}#end f_ReadData<-function(d)
###############################################
f_output_minisum<-function(ls_result_minimsum,s_output_file ){

v_Call_mkt_1=unlist(ls_result_minimsum[1])
v_Call_mkt_2=unlist(ls_result_minimsum[2])
v_maturity_dt=unlist(ls_result_minimsum[3])
v_X1=unlist(ls_result_minimsum[4])
v_X2=unlist(ls_result_minimsum[5])
v_S=unlist(ls_result_minimsum[6])
v_risk_free=unlist(ls_result_minimsum[7])
v_dividend=unlist(ls_result_minimsum[8])

v_Minisum=unlist(ls_result_minimsum[9])
v_sigma_star=unlist(ls_result_minimsum[10])
v_flag=unlist(ls_result_minimsum[11])
v_fi=unlist(ls_result_minimsum[12])


m_result_minimsum=rbind(v_Call_mkt_1,v_Call_mkt_2,v_maturity_dt,v_X1,v_X2,v_S,v_risk_free,v_dividend,v_Minisum,v_sigma_star,v_flag, v_fi )


write.table(m_result_minimsum , file = s_output_file, sep = ",", col.names = NA, qmethod = "double")
print(m_result_minimsum )


return(0)
}# end f_output(ls_result_minimsum,s_output_file )


#############################################
###I.2. MAIN PART of the PROGRAM MINISUM
#############################################
# read data
s_input_file="C:/Users/diepnguyen/Documents/MSF- Clark/Implied Binomial Tree//Research/R/SampleData_SK_012011.csv"
ls_data=f_ReadData(s_input_file)
print(ls_data)

#compute minisum

  ls_result_minimsum=f_minisum(ls_data)
print(ls_result_minimsum)

#print out outcome
   s_output_file = "Result_minisum.csv"
   f_output_minisum(ls_result_minimsum,s_output_file )




getwd() #See directory to the output file
